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Abstract 

Dirichlet Process Mixture (DPM) models have been increasingly employed to 
specify random partition models that take into account possible patterns within 
the covariates. Furthermore, to deal with large numbers of covariates, methods 
for selecting the most important covariates have been proposed. Commonly, the 
covariates are chosen either for their importance in determining the clustering of the 
observations or for their effect on the level of a response variable (when a regression 
model is specified). Typically both strategies involve the specification of latent 
indicators that regulate the inclusion of the covariates in the model. Common 
examples involve the use of spike and slab prior distributions. In this work we 
review the most relevant DPM models that include covariate information in the 
induced partition of the observations and we focus on available variable selection 
techniques for these models. We highlight the main features of each model and 
demonstrate them in simulations and in a real data application. 

Keywords. Dirichlet Process Mixture models, random partition models, Bayesian 
variable selection, spike and slab distributions, model misspecification. 


1 Introduction 

Bayesian nonparametric literature has been increasingly focusing on models that can 
cluster observed units according to possible patterns in the covariate space. A common 
strategy is usually referred to as Random Partition Model with Covariates (RPMx, Muller 
and Quintana [2010]) and has been successfully applied to a wide range of real-data prob¬ 
lems, including epidemiology (Park and Dunson [2010]), survival analysis (Muller et al. 
[2011]), genomics (Papathonias et al. [2012]), pharmacokinetics and pharmacodynamics 
(Muller and Rosner [1997]), finance (Griffin and Steel [2006]). 

Usually, an RPMx is constructed starting with a Dirichlet Process Mixture (DPM, 
Lo [1984]) model. This is characterized by specifying a Dirichlet Process (DP, Ferguson 
[1973]) prior on the parameters of the sampling model. The popularity of these models 
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is due to fact that they allow for high flexibility and that the posterior distribution of 
interest can be explored by efficient computational algorithms. DPM models induce a par¬ 
tition of the observations in clusters, with the probability of belonging to a specific cluster 
proportional to the cluster’s cardinality a priori. This imposes a normal behavior on the 
partition. Recently, a wealth of research has been focussing on enriching the clustering 
structure, by introducing dependence of the cluster probability on covariates. Moreover, 
RPMx have been extended to embed latent parameters with the aim of performing vari¬ 
able selection. The role of the latent variables in the RPMx framework consists primarily 
in identifying the subset of variables that are more discriminant in terms of the partition. 
The variable selection output is just the posterior distribution of the latent indicators, 
which are commonly treated as any other model parameter and whose distribution is 
often approximated by Markov Chain Monte Carlo (MCMC) techniques. 

The main objective of this paper is to review the most relevant RPMx models, defined 
through Dirichlet Process Mixtures. We dedicate particular attention to the available 
variable selection techniques. The rest of the work is organized as follows. In Section 2 
we review the relevant theory about DPM models. In Section 3 we present the relevant 
literature about DPM with covariates, while in Section 4 we review available variable 
selection methods. In Section 5 we present a simulation study and in Section 6 results of 
a real application are shown. We conclude with a final discussion in Section 7. 


2 Dirichlet Process Mixture Models 


The Dirichlet Process (DP) is a distribution over random distributions (Antoniak [1974], 
Ferguson [1973]). A constructive definition is presented by Sethuraman [1994], who 
showed that if a random probability measure G is distributed according to a DP with 
precision a G M"*" and center measure Go defined on the metric space 0, then 


G - '^^pkSek, 

k=l 


( 1 ) 


where the elements 61 , 62 , ■ ■ ■ are iid realizations from Go, is the Dirac measure that 
assigns a unitary mass of probability in correspondence of location 6 k and the V’fc cire 
constructed following the stick breaking procedure (see Ishwaran and James [2001] for 
details); 

k-l 

'ipk = (l)kW0-- ( 2 ) 


with (pk ~ Beta(l,Q;). By construction 0 < V’fc < 1 and V’fc ~ 1- The resulting 

random probability measure G is defined on the same support of Go, i-e. 0. A more 
compact notation is G ~ DP (a. Go). 

Another common representation of the DP, which allows for efficient MCMC schemes, 
has been provided by Blackwell and MacQueen [1973]. Let us consider a sample of n 
components 0 = [61,... ,6^) from a random distribution G. If G is distributed as a 
DP(a, Go), then by integrating out G from the joint distribution of ^i, ..., 6 n, we obtain 
the predictive prior distribution of 6 i given which is the vector obtained by removing 
the i-th component from 0 \ 


6 , I 


a + n 


3T + 


a 




a -f n — 1 


-Go. 


(3) 
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Equation (3) is generally referred to as the Blackwell-MacQueen urn scheme. In 
particular, the first component of Equation (3) can be rewritten as + 

n — 1), where rij is the number of observations that have value equal to 6 *. The vector 

g{i)* ^ 

contains the unique values of the sequence Since Equation 
(3) is a mixture of atoms and of a diffuse measure, there is a positive probability that 
k < (n — 1). This aspect is due to the discreetness of the DP samples (Blackwell [1973]): 
there is a positive probability of ties, i.e. that two random draws from G ~ DP(-, •) are 
identical. From Equation (3) it is also clear that there is a higher probability that a new 
(as yet unobserved) unit will be assigned to a larger cluster (in terms of cardinality). 

This aggregating property of DP makes it particularly effective to deal with clustering 
problems. In fact, arguably the most famous application of the DP is the Dirichlet Process 
Mixture (DPM) model (Escobar and West [1995], Lo [1984]), a class of models that can 
be expressed hierarchically as follows: 

yi,...,yn\0i,...,9n p{yi I 9i) 

9u...,9n\G ~ G (4) 

G ~ DP(a,Go). 

This model assumes individual level parameters 0,, for i = l,...,n. Throughout the 
paper we use the word model to indicate the joint probability distribution of all unknowns, 
including data and parameters. With a slight abuse of terminology we use model and 
method interchangeably. The vector of parameters will have some ties with probability 
greater than zero. This is because we set each one of them to have a distribution G 
which is a DP. This will have two main consequences: (i) the sequence 6 = ( 6 ^ 1 , ..., 9n) 
reduces to the sequence of its unique values 6 * = (6*]',..., 0^), with k < n, (ii) the vector 
s = {si,, Sn) with Si G {l,...,/c}, which associates each observation with a specific 
value among the components of the vector 6 *, defines a partition of the observations. In 
practice, the sets of this partition can be interpreted as clusters of individuals. 

An alternative representation of the DPM model is given by: 

yi,...,yn\G ~ p{y \ G) 

p(y I G) = fp{y\ e)G(de) (5) 

G ~ DP(a,Go). 

Recalling the discrete nature of the DP samples as well as its representation in Equation 
( 1 ), we can rewrite the sampling model as an infinite mixture model: 

OO 

2 / 1 , • • •, 2 /n I G ~ ^ 'ipkPiy I 9k). 

k=l 

Let pn denote the partition of the n observations implied by s. It is easy to prove 
that the prior distribution for induced by the DP prior is: 

= ( 6 ) 

where = a{a + 1)... (a + n — 1) (take i in Equation (3) to be the last observation 
for i = 1,... ,n, exploiting the exchangeability of the Blackwell-MacQueen urn). This 
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defines an Exchangeable Partition Probability Function (EPPF, see Pitman [1996]), where 
exchangeability arises from the fact that the partition does not depend on the labels of 
the observations or of the clusters, but only on the cardinality of the groups. 

Therefore, a DPM model can be represented as a Random Partition Model (RPM, 
see Lau and Green [2007] for details) through p{pn)- An RPM is characterized by within- 
cluster-submodels and by a prior distribution on the partition. This is evident when 
writing the joint probability model of the DPM in Equation (4) (Lo [1984]) as : 

p{pn,y,6*) oc JJ < JJ [pivi I e*)] go{6*)a{nj - 1)! i , (7) 

j=l [ i&Sj J 

where go is the density associated with the distribution Go, while Sj = {i : Si = j, for i = 
1,... ,n}. Compared with Equation (4), this is the joint density with G integrated out 
and reparameterized in terms of the partition and unique values. The term a{nj — 1)! is 
called cohesion function for the j-th group and is denoted by c{Sj). Since p{pn) can be 
seen as the product of the cohesion functions for each of the groups, this links the DPM 
with a specihc type of RPM called Product Partition Model (PPM, Barry and Hartigan 
[1992], Hartigan [1990]), characterized in the same way. 

Extensions to the model in Equation 4 and 7 can be achieved by employing more 
general classes of prior distributions for G. For a detailed review see Lijoi and Priinster 
[ 2010 ]. 

Using Equation (3), it is possible to specify the conditional posterior distribution of 
6i for the model in Equation (7) as follows: 

p{9i I 6^^\y) oc 'Y^piVi I 9i)5eX9i) + «/ viVi I 9)dGo{9)go{9i \ yi). (8) 

Particularly within a regression framework, recent Bayesian literature has focussed 
on dehning RPM allowing for covariate information when inferring the partition of the 
observations. This can be obtained by modifying the cohesion function to account for 
covariates patterns. At the same time, there has been an increasing interest in perform¬ 
ing variable selection within the context of RPM with covariates to identify the most 
informative variables for the partition. In the next sections we will first review the main 
methodologies for specifying DPM-based RPM with covariates and then we will present 
state of the art procedures for variable selection. 

3 Covariate dependent DPM 

Let us consider a matrix of covariates X with n rows and D columns and let Xi denote 
the th row. In many applications, it is desirable to express a prior distribution on the 
partition that is a function of X, i.e. p{pn \ X), instead of letting the probability of 
the partition depending (a priori) only on the cardinality of the clusters. This type of 
models has been called Random Partition Model with Covariates (RPMx). See Muller 
and Quintana [2010] and Dunson [2010] for surveys. 

We will focus on RPMx that admit a product partition representation {e.g. the DPM 
models). To this end we allow the cohesion function to include the covariates, however, 
we assume that the overall probability of a partition is still specihed as the product of 
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cohesion functions of each cluster: 

k 

p(Pn\X)a:'[{c(Sj,X'’-), ( 9 ) 

where, for j = 1,... fc, Xj" is the subset of the rows of X associated with cluster j of 
the partition In the following sections we will review the most popular choices of 
covariate dependent cohesion functions. 

In a regression framework, when the research interest is in modeling the relationship 
between a response variable y and a set of covariates X, i.e. studying the density of 
p{y\x,o), the application of RPMx models has been very frequent. This is mainly for 
two reasons. First, RPMx are flexible models, which allow to cluster the observations 
according to patterns within the covariates and then to specify a cluster-specihc regression 
model. Secondly, they often lead to improved predictions: if we want to predict the 
response for a new subject with a specihc set of covariates, then a RPMx model will 
assign higher probability that the new subject belongs to the cluster that contains the 
most similar covariate prohles. 

3.1 Augmented Response Models 

The most common strategy to include information about X into the partition model in a 
DPM framework has been to treat each covariate as a random variable, i. e. by specifying a 
suitable probability model. Muller et ah [1996] were the hrst to introduce this idea within 
the DPM framework. In their work they consider an augmented model dehned on Z = 
[y.x) and their objective is to estimate the smooth function g{X) = E{y \ X). They 
approach the problem by modeling as a DPM of {R + D)-dimensional distributions, 
where R is the dimension of the response variable (usually R = 1). Let A* be the matrix 
containing the unique parameters for the k clusters, (A]|,..., A|,). Considering now a new 
observation z = {y,x), its predictive distribution can be derived as: 

k 

Piy, * I M) oc ^ njp{y, x \ A*) + a 

Assuming uncertainty about the realized value of x, which might be a reasonable and 
necessary assumption when x is measured with error or not exactly known in real appli¬ 
cations, allows us to rearrange the latter equation as 

k 

P{y I A*) oc '^njp{x I A*)p{y \ x,A*) + a 
i=i 

using Bayes’ theorem. The quantity njp{x \ A*) depends on the cardinality of group j 
and on a measure of how likely it is that the new observation will be clustered in group 
j, based on the value of its covariates. The latter is the likelihood of the observed x. 
The smooth function g{X) is then estimated by taking the expectation with respect to 
p{y I x,A*). Muller, Erkanli and West describe in details the case where Z is a mixture 
of multivariate Gaussian distribution, which leads to simplihed calculations for g{X). 

A similar approach has been adopted by Miiller et ah [2011]. They originally propose 
a modihcation of a PPM, the PPMx (PPM with covariates), to incorporate measures of 


jp{y I x,A)p{x I A)dGo(A), 


p{y,x I A)(iGo(A). 
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similarity among the covariates within each cluster employing the following structure for 
the prior of the partition of the observations: 

k 

P(pj (10) 

J = 1 

where /(•), called similarity function, is an ad hoc function that takes large values for 
highly similar values of the covariates. The authors propose as a default choice to specify 
/(■) as a probability density. They show under mild conditions that f{Xj") can be 
seen as the likelihood of the covariates belonging to cluster j, from which the cluster 
specihc parameters have been integrated out. Given the cluster specihc parameters for 
the covariates, the joint probability of a PPMx is: 

f{y,X,0*,Cl,...,CD,Pn)^ 

k 

n n I I Cl, 

j=l ieSj 

where 6* and Ci, • • •, Cd include the unique values of the parameters of the distribution 
of the response and of the covariates for the k clusters respectively. Equation (11) shows 
that the PPMx is a generalization of the methodology proposed in Muller et al. [1996]. 
Taking c{Sj) in Equation (11) to be the cohesion function implied by the DP and the 
covariates to be random variables with distribution p{xi \ Qi,... ,( 10 ) (thus allowing 
the similarity function to be a valid probability density for the covariates), the PPMx 
simply reduces to a DPM on the joint distribution of the response and the covariates 
representable by the following hierarchy: 



Vly ■ ■ ■ 1 Vn 1 

x,e 

ind 

p{yi 1 Xi,9i) 


Xi, . . 

• , 1 Cl, ■ ' 

• - , Cn 

ind 

P{Xi 1 Ci) 

(12) 

XiXi, 

),•••, {9ni Cr 


iid 

rsj 

G 




G 

rsj 

DP(q:, Go), 



with Go = Goe x Gqc, O = and Ci = (0i,---,0d)- Both Equation (11) 

and (12) dehne a PPMx, in which 6 and are assumed a priori locally independent 
but globally dependent. Therefore, every DPM can be represented as a PPMx, but the 
reverse is not always true. For this relation to hold, it is necessary that p{yi, Xi \ 9i, Ci) = 
p{yi I 9i,Xi)p{xi I Ci)- III fliis perspective the PPMx generalizes the work by Miiller 
et al. [1996] allowing for the possibility of user-specihc models for the covariates (via the 
similarity function). An example of PPMx is represented by the work of Barcella et al. 
[2015] which specihes a model for dealing with binary covariates containing information 
about symptom prohles. 

Alternatively, Park and Dunson [2010] have proposed the Generalized Product Par¬ 
tition Model (GPPM). The authors discuss how to incorporate covariate information in 
the conditional prior distribution in Equation 3. This results in a generalized Polya urn 
scheme from which they derive a covariate dependent version of the PPM which show 
the same joint model in Equation (11). 

Within the PPMx framework in Equation (11), the sampling model p{yi \ 9*,Xi) 
does not necessarily need to be a linear regression. Hannah et al. [2011] have extended 




( 11 ) 
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Equation (11) to the broader Generalized Linear Model (GLM) framework through the 
appropriate specification of p{yi \ 6*,Xi). This generalization allows the user to handle 
different types of data. They refer to this model as DP-GLM (see also Shahbaba and 
Neal [2009]). A parametric version, i.e. with a finite number of mixture components, of 
the DP-GLM constitutes a particular case of the Hierarchical Mixture of Experts (HME) 
model introduced by Jordan and Jacobs [1994] and specified in a Bayesian framework by 
Bishop and Svenskn [2002]. 

Profile Regression (PR; Molitor et ah [2010]) is another prominent example of aug¬ 
mented response models. In the original formulation this model handles a binary outcome 
y = (i/i,... ^Un) which is common in epidemiological applications, but the model is eas¬ 
ily generalized to different types of response variable. The PR model consists of two 
submodels. The first one is the model for the response: 

l/i I Pi ~ Bernoulli (pi), 
with a logistic regression on the mean pp 

=^i + (13) 

where Wi is a set of confounding variables with coefficients k while 6i is an individual 
random intercept. 

The second submodel is a mixture model on the covariates, such that conditioning on 
the cluster assignment vector, the probability of a specific covariate profile becomes: 

Xi\ Ci ^ p{xi\ Ci)- (14) 

When Xi is a vector with D components, we can model each component independently 
and we can treat Q as a vector containing the parameters for each component of the 
profile, i.e. Ci = (Oi, • • •, Od)- 

In order to consistently estimate the posterior distribution of the partition and the 
partition specific parameters, the authors propose to model jointly the random intercepts 
in Equation (13) and the parameters of the covariates sub-model in Equation (14) ac¬ 
cording to an unknown distribution G, which follows a DP with parameter a and Go, 
with Go being the product measure of Goe and Go(^. Expressing the joint model in terms 
of the implied partition and the cluster-specific parameters, the PR can be equivalently 
represented as the PPMx in Equation (11). 

For the augmented response class of models, R packages are available for the PPMx 
(https://www.ma.utexas.edU/users/pmueller/prog.html#PPMx) and for PR (http: 
//cran.r-project.org/web/packages/PReMiuM/). 

3.2 Dependent Dirichlet Process 

An alternative way to include covariate information in DPM is to allow the weights and/or 
the locations in the stick breaking construction of the DP in Equation (1) to depend on 
covariates. In particular, a this can be represented in the following way: 

OO 

Gx = '^'ipk{x)Se^ (15) 

k=l 

k-1 

ijkix) = (j)k{x)Y[[i - 4>jix)], 

1=1 


7 



under the constraint that = 1- V’fc(') is a function of the covariates. In this 

context X represents a point in some covariate space X and 4>k{x) is a realization of a 
Beta distribution with parameters equal to 1 and a{x), the latter being the (positive) 
realization of a stochastic process indexed at x E X. The model dehned in Equation 
(15) is a particular case of the Dependent Dirichlet Process (DDP, MacEachern [1999]). 
Each Gx is still marginally a DP for each x. In its original formulation, the DDP model 
allows for both covariates dependent weights (as in Equation (15)) as well as for covariates 
dependent locations. In many applications the original formulation of the DDP has been 
reduced to accommodate covariate dependent locations only (examples are De lorio et al. 
[2009, 2004], Gelfand et al. [2005], Duan et al. [2007], among others). However in terms 
of the random partition models. Equation (15) (or the version including additionally 
covariate dependent locations) presents the most relevant construction (see Muller and 
Quintana [2010]). In this case the specihcation of a distribution for 'ipk{x) is central, as it 
determines the structure of the dependence between the covariates and the weights, and 
consequently the way the covariate prohles inform the clustering structure. 

Although assuming 'ipki^) sue product of Beta random variables guarantees that the 
Gx are marginally (for each level of the covariates) a Dirichlet process or other known 
processes (see, for example, the covariates order-based stick-breaking of Griffin and Steel 
[2006]), several authors have preferred to use different models for (pki^) in order to allow 
for more flexible stick-breaking processes. The resulting processes do not belong to DDP 
anymore. Examples include the kernel stick-breaking {i.e. when 0(-) is user-dehned 
function with codomain in (0,1) which often captures the distance of the covariates from 
centroids) in Reich and Fuentes [2007], Ghnng and Dunson [2011], Dunson and Park 
[2008] and Griffin and Steel [2010], the probit stick-breaking [i.e. 0(-) is the cumulative 
distribution function of a Normal density, whose input can be a function of the covariates 
or alternatively a spatial process indexed to the covariates) in Rodriguez et al. [2009], 
Ghnng and Dunson [2009], Rodriguez and Dunson [2011] and Arbel et al. [2016] and the 
logistic stick-breaking [i.e. is a logit function, whose argument is a function of the 
covariates) in Ren et al. [2011] among others. See Foti and Williamson [2015] for a review. 

The choice of the distribution for 'ipk{x) determines the DDP (or a dependent stick¬ 
breaking), which can then be used as mixing measure in a hierarchical model leading 
to; 

OO 

PiVi \x,6) = ^V’fc(aj)p(2/i I Ok)- 

k=l 

Note that it is also possible to assume a regression sampling model for y, i.e. p{yi \ x, 6k) 
instead of p{yi \ 9k). 

A related approach is the Weighed Mixture of DP (WMDP) by Dunson et al. [2007], 
which can be thought of as a hnite mixture of DP distributed components, one for each 
covariate level. The weights of this mixture are specihed as functions of the covariates. 
The resulting random measures maintain covariate independent locations and can be used 
conveniently to specify an inhnite mixture model with covariates dependent weights. 

3.3 Other Methods 

In this section we briefly present two other methods that can be used to specify covariate 
dependent DPM. 

The hrst one is the Restricted DPM (RDPM) model introduced by Wade et al. [2013]. 
The authors modify the usual structure of the DPM models by imposing restrictions to the 
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distribution of the partition of the observations to follow the covariate proximity. For ex¬ 
ample, let us consider n instances of a univariate covariate, xi,..., and the permutation 
of 1 ,..., n given by ordering increasingly the covariate values, namely cr 2 ;(l),..., (Tx{n). 
The RDPM restricts the prior probability over the partition of the observations implied 
by a DPM and considers only the partitions for which < So-^(n). It can 

be shown that this construction satishes the Ewens sampling law (Ewens [1972]) for the 
probability of the cluster frequencies. This same law is satished by partitions implied by 
Equation (3). This class of models is appealing because it does not assume any distri¬ 
bution on the covariates when accounting for the covariate similarity. The authors show 
how to perform posterior inference in the RDPM through efficient MCMC algorithms. 
The mixing properties of the MCMC scheme are improved by restricting the support of 
the random partition. 

A second alternative is represented by the Enriched Dirichlet Process Mixture (EDPM) 
model described in Wade et ah [2014]. The strength of this method consists in its abil¬ 
ity to create nested partitions (he. partitions within sets of a partition). To this end, 
the authors specify a DPM model for the response variable, setting a DP prior on the 
parameters of the sampling model for y. A DP prior, dependent on the parameters of 
the response, is used for the parameters of the sampling model on the covariates. This 
construction leads to a nested clustering structure of the observations; a hrst level of 
clustering is at the response level, whereas a second level is obtained within the clusters 
formed at the hrst stage according to a DPM model on the covariates. 

3.4 Remarks 

Covariate dependent Dirichlet process mixture models have been increasingly used in 
practice, especially when the objective is to specify hexible regression models. The main 
motivation underlying the use of such models is to improve predictions, in comparison to 
other possible nonparametric cluster-wise regression models. The latter has been demon¬ 
strated in simulation for augmented response models in Cruz-Marcelo et al. [2013]. The 
improvement in predictions is the result of substituting the traditional mixture weights 
in DPM models, which depend on the cardinalities of each cluster, with some function of 
the covariates. In this way the relation between covariates and response is studied within 
clusters of observations, whose assignment probabilities vary across the covariate space. 

The review of covariate dependent Dirichlet processes presented in this section shows 
that there are mainly two strategies for specifying such models in the context of Dirichlet 
processes. The hrst way consists of modeling jointly the response and the covariates as 
a Dirichlet process mixture of multivariate distributions. The main advantage of using 
this technique is its computational simplicity. In fact, for all types of covariates the main 
model remains a DPM, which has computational advantages allowing the use of efficient 
algorithms by MacEachern and Muller [1998], Neal [2000] for posterior inference. For 
these models it is also possible to integrate out the variability on the mixing measure 
so that the conditional prior distributions on the parameters of the mixture model can 
be expressed as a modihed Blackwell-MacQueen urn which includes the covariates (see 
Park and Dunson [2010]). On the other hand, the main disadvantage of this strategy 
is related to the fact that for high dimensional covariate space the likelihood of the 
augmented response variables becomes dominated by the portion relative to the covariates 
and consequently the response does not inform effectively the clustering. 

The second technique relies on modifying the stick-breaking process through which 
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the weights of the traditional DPM models are constructed to include covariates. All 
contributions to this held can be divided between those that assume DPM models for 
each level of the covariates and those which do not. In the hrst case the stick-breaking 
procedure at each covariate level has to involve a sequence of Beta(l, a) random vari¬ 
ables. This may be a limitation in incorporating complicated covariate dependences in 
the weights, thus stick-breaking procedures which involve link functions that map some 
regression of the covariates into the (0,1) set have progressively been employed. Once 
a convenient link function is found, a variety of types of dependence can be accommo¬ 
dated in the weights, which is the main advantage of these techniques. However, this 
kind of models often leads to poor inference when few observations are available for each 
covariate level (even more so in presence of continuous covariates). Furthermore, poste¬ 
rior inference may require more sophisticated algorithms (as the slice sampler by Walker 
[2007] or retrospective sampler by Papaspiliopoulos and Roberts [2008]) or truncation of 
the inhnite mixture to some fixed level for allowing the use of the blocked Gibbs sampler 
by Ishwaran and Zarepour [2000]. 

4 Covariate Dependent DPM and Variable Selection 

Increasing research interest has been devoted to develop variable selection strategies in 
covariate dependent DPM models. Bayesian methods for variable selection have a long 
history and a variety of different techniques have been proposed to achieve this task (see 
O’Hara et al. [2009]). Within the regression framework, this corresponds to evaluate the 
uncertainty about the selection of covariates to include in the model. One of the most 
common way to perform Bayesian variable selection in regression framework consists in 
specifying prior distributions favoring shrinkage toward zero on the regression coefficients. 
Similarly, indicators can be included in the model to select which covariates are active in 
the model. Alternatively, a prior distribution directly over the model structure can be 
specihed. In this section we describe exclusively variable selection techniques proposed for 
covariate dependent DPM models. We deal separately with tools for augmented response 
models and Dependent Dirichlet Process. 

4.1 Variable Selection for Augmented Response Models 

Product Partition Model with Covariates (PPMx) 

A variable selection strategy for the PPMx has been proposed by Muller et al. [2011] 
and described in details by Quintana et al. [2015]. Without loss of generality we start 
our discussion by considering the PPMx from the RPM point of view. R is possible to 
rewrite the similarity function in Equation (10) as the product of the similarity functions 
of each individual covariate, i.e. = Y[d=ifi^jd)^ where Xj2 is the sub-vector 

of elements of column d of X which includes the elements corresponding to cluster j. 
Variable selection is then introduced employing binary indicators 7 *^ for j = 1,... k and 
d = 1,..., D within the distribution of the partition: 

k D 

p{pn I X,7) c^Y[c{Sj)Y[f{xP2yU. (16) 

i=l d=i 

The presence of the binary indicators allows the probability of the partition to depend 
on a subset of covariates within each cluster. In fact, 7 *^^ = 0 eliminates the effect on 
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the distribution of the partition of covariate d in cluster j. In this setting, extra care is 
required for the specihcation of /(■). In order to perform variable selection, /(■) must 
always take values larger than 1 (otherwise excluding a covariate always increases the prior 
probability). The authors discuss convenient choices of /(■). The model is completed by 
introducing in the hierarchy a prior distribution for the indicators. In particular, the 
authors propose to use a Bernoulli prior distribution assuming a logistic link for the 
probability of success. 

A different method for performing variable selection in PPMx framework is presented 
by Barcella et ah [2015], which extends the work of Kim et al. [2009] to the augmented 
response class of models. The authors specify a joint DP prior on the regression coeffi¬ 
cients and the parameters governing the distribution of the covariates, assuming a priori 
local independence between the two sets of parameters. Assuming a spike and slab base 
measure for the regression coefficients, this model allows to perform cluster specihc vari¬ 
able selection, while, at the same time, the clustering structure is informed by both the 
covariate prohles and the relationship between response and covariates. The authors refer 
to this model as Random Partition Model with covariate Selection (RPMS). 

More formally, the RPMS can be represented by a hierarchy similar to the one in 
Equation (12): 

yi,... ,yn \ X,&,X Normal(?/i I , A) 

n D 

nn Bernoulli(xjd | Qd) (17) 

i=l d=l 

(e,.C,),...,(e„,C)|G " G 

G ~ DP(a,G'o), 

where © and Z are matrices of parameters with n rows and D columns. For i = 1,..., n, 
/3j is a D-dimensional vector and is a row of ©; similarly, is a D-dimensional vector 
and a row of Z. The RPMS in Equation (17) is designed in the original formulation to 
handle binary covariates, even though changing the specihcation of the distribution of 
the covariates enables us to include different types of variables. The center measure Gq 
has the following form; 

D 

Go = JJ{[7rA(6'd) + (1 - '^d)N{9d \ /id, rd)]Beta(Cd | &c)}> (^8) 

d=l 

and we can rewrite Go = Goe x Go(^- Following Kim et al. [2009], Barcella et al. in¬ 
duce super-sparsity to the matrix of the regression coefficients following the hyperpriors 
structure presented by Lucas et al. [2006]. 

Additionally, in PPMx framework Kunihama and Dunson [2014] consider an aug¬ 
mented response model and they propose a method for testing for conditional indepen¬ 
dence of the response and a specihc covariate given all the other covariates. This involves 
the conditional mutual information for measuring the intensity of the dependence. 

Profile Regression (PR) 

Papathomas et al. [2012] investigate the problem of performing variable selection within 
the Prohle Regression framework when all the covariates are categorical (see also Pa¬ 
pathomas and Richardson [2014]). Let us recall that PR can be decomposed into two 
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sub-models; a model on the covariates and one on the response. These are linked by 
using a joint DP prior on the set of parameters common to both the submodels. In order 
to introduce variable selection we need to rewrite Equation (14) in the following way: 

D 

d=l 

Variable selection is then performed by replacing the distribution of each covariate with: 

P^^iXid I Cjd, TTd) = 7ldP{Xid I Cjd) + (1 - T^d)rdiXid), (19) 

where the superscript VS indicates that the implied probability has been modified to 
perform variable selection, G ( 0 , 1 ) is a continuous weight and rd{xid) indicates the 
proportion of times covariate d takes value Xid- From Equation (19) it is evident that 
large values of vr^ indicate that covariate d is informative in terms of clustering. In this 
setting a Beta hyperprior distribution for each tt^ or alternatively a mixture of a Beta 
distribution and Dirac measure (with Bernoulli distributed indicators) may be preferred 
to induce extra sparsity. The authors compared their approach that uses continuous 
weights to a version that employs cluster specific binary indicators for each covariate. 
The latter idea can be represented in the following way: 

P^^^ixid I C]dn*d) =p{xid I C]dy^‘^rd{xidY^~'^^'^\ 

where 7 *^ = 1 indicates that covariate d is informative with respect to cluster j. This 
approach is a generalization to Profile Regression of a solution proposed by Chung and 
Dunson [2009]. In contrast with the continuous case, the natural choice of prior dis¬ 
tribution for each 7 *^ is Bernoulli with mean distributed as a Beta distribution. Extra 
sparsity can be achieved substituting the latter Beta distribution with a mixture of a 
Beta distribution and Dirac measure (with Bernoulli distributed indicators). 

The results presented by Papathomas et al. [2012] and obtained employing the extra 
sparsity alternative of both variable selection methods described above show comparable 
performances of the two methods in terms of variable selection, although preference is 
given to continuous weights due to faster MCMC convergence. 

An extension of the methods above has been proposed by Liverani et al. [2015] to 
deal with continuous covariates. This consists in modifying Equation (19) substituting 
T^dixid) with a suitable summary, for example the observed mean of the d-th covariate. 

4.2 Variable Selection for DDP 

To the best our knowledge, general variable selection strategies have not been imple¬ 
mented in the DDP framework. However, in the case of the dependent stick-breaking 
process Chung and Dunson [2009] show how to perform covariate selection when the 
weights of the random probability measure are constructed by a probit link stick-breaking. 
Recalling the stick-breaking procedure in Equation (15) the following specification is pro¬ 
posed: 

00 

Gx = '^'ipk{x)6e^ ( 20 ) 

k=l 

k-1 

'ipki^) = 4>(z/fc(a;)) JJ[l-<h(i/j(a;))], 

1=1 
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where $(■) is the standard normal distribntion and is a predictor which can be 
specified for example as Vk{'^) = Variable selection is then achieved by introdncing 
binary indicators: 

D 

ik ~ !!?(«“ I (21) 

d=l 

where ad denotes the covariate specific parameters of the distribntions of ^kd for all k. 
Considering a regression sampling model p{yi \ Xi, 6k), it is possible to link the resnlts of 
the variable selection performed in Eqnation (21) directly to the parameters 9kd fo fhe 
regression model for the response so that when jkd = 0 both 9kd and ^kd are set eqnal to 
0 . 

4.3 Remarks 

In this section we have reviewed the available methodologies for performing variable 
selection in covariate dependent random partition models. We conld distingnish between 
two main approaches: one selects the covariates for their importance in terms of clnstering 
{e.g. the variable selection methods proposed for the PPMx or for PR) and one selects the 
covariates which are relevant for explaining the level of the response within each clnster 
when a regression is specified for the model of the response {e.g. RPMS). 

When a regression sampling model is employed none of the approaches above allows 
in principle to exclnde a covariate from the model. For example, if the RPMS exclndes a 
covariate as infinential on the level of the response in a certain cluster, however it cannot 
exclude the same covariate from affecting the clustering. Similarly, in PPMx excluding 
a covariate from affecting the clustering does not imply automatic exclusion of the same 
covariate from the regression sampling model. 

A more elaborate solution which links variable selection in terms of clustering and 
association with the response level has been presented by Chung and Dunson [2009]. 
This proposal employs common binary indicators for each covariate in both the sampling 
model and the model of the weights. This implies that if a covariate is excluded from the 
model of the weights is automatically excluded from the model of the response. 

5 Simulation Study 

Specifying covariate dependent weights in mixture models has the advantage of making 
posterior inference robust to model misspecification. We illustrate this point with two 
simulation studies (one presented in Supplementary Material) in which we compare the 
results of the Random Partition Model with covariate Selection (RPMS, Barcella et al. 
[2015]), the Profile Regression (PR, Molitor et al. [2010]) the Probit Stick Breaking Pro¬ 
cess Mixture Model (PSBP-MM, Chung and Dunson [2009]) and the model described in 
Kim et al. [2009], which, for simplicity, we refer to as Spike and Slab Model (SSM). The 
latter model is simply a DPM model of regressions for which the center measure of the 
DP is chosen to be a spike and slab distribution similar to the one adopted in the RPMS. 
In other words, RPMS and SSM have the same hierarchical structure of Equation (17), 
except for the model on the covariates. We use for both the RPMS and SSM the same 
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hyperpriors for the center measure of the DP: 


TTi, . . . , 71 d I 1 ^ 1 , ■ ■ ■ , 

(Ul, . . . , Ud 

Ti,. . . ,Td 


D 

~ Ilid - UdjSoiTid) + a;dBeta(7rd | a^, b^)) 

d=i 

D 

~ Y\_^eta{ud\ auj,b^) 
d=i 
D 

n Gamma(rd | ar,br). 

d=i 


( 22 ) 


For PSBP-MM we specify the following sampling model: 


Vil Gx ^ 


j NoTmal{x0'^, X)dGx{d), 


where Gx is the process described in Equation (20) and Equation (21), for which we 
assume Uki^) = and 


D 

^k\jk ~ JjNormal(^M I 

d=l 

jkd I i^d ~ Bernoulli(7fcd | Kd) 

Kd\ud ~ Beta(Kd I a«,6«)“‘'((5o(fi;d))^^"“''^ 
Ud ~ Bernoulli(0.5). 


We mostly follow the specihcation of the PSBP-MM given in Chung and Dunson [2009], 
but without modeling the distance of the covariate values from the centroids within the 
weights. Furthermore, the center measure of Gx{0) is assumed to be the product of 
D independent distributions: Normal distributions with mean 0 and precision Td for 
the included covariates and Dirac measures located at 0 for the non-included covariates. 
We employ Gamma prior distributions with parameters a,- and br for each ti ... ,td ■ 
Chung and Dunson [2009] use instead a multivariate Normal distribution for the included 
covariates, assuming mixtures of gf-priors (Liang et al. [2008]) on the covariance matrix. 

In what follows, we only consider binary covariates. We choose the same hyperpa¬ 
rameters for the RPMS and the SSM models: 0 ,^ = = 0.15, a^j = l,b^ = 0.15,0,- = 

br = I, a\ = b\ = lytta = ba = 1 and, only for the RPMS, = 6 ^ = 1. We do not 
update the parameter fid and we hx it equal to 0 for all d. As mentioned above, in both 
cases posterior inference is performed through MCMC algorithms. We initialize the al¬ 
gorithm starting with one cluster and hxing the regression coefficients equal to zero and 
the parameters for the covariates equal to 0.5 (this last specihcation is required only for 
the RPMS). We run 15000 iterations, discarding the hrst 5000 as burn in. 

The PR has been initialized with the default values of the R package PReMiuM and 
10000 samples have been saved after discarding the hrst 5000. A Normal distribution 
for the response and a Bernoulli distribution for the covariates have been assumed. Con¬ 
founding variables have been ignored. We focus exclusively on the variable selection via 
continuous indicators (see Equation (19)) following the suggestion of the authors. 

Finally, the hyperparameters for the PSBP-MM have been hxed to the following 
values: 0^ = 5^ = 0.5, ar = 1 and br = 5, = 0 and = 0.1. Following Chung and 

Dunson [2009] , posterior inference has been performed using a blocked Gibbs sampler (see 
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Ishwaran and James [2001]) which requires a truncation level K for the inhnite mixture 
model. We £x = 20. We run 15000 iterations, discarding the hrst 5000 as burn in. 

Convergence of the chains have been investigated by trace plots and computing au¬ 
tocorrelations of the samples. The results show evidence of convergence for the chains of 
all estimated parameters. 

5.1 Cluster-wise linear regression model with interactions 

We simulated a dataset with n = 200 observations. We considered two binary covariates 
{i.e. D = 2) and each entry Xid of the design matrix was generated from a Bernoulli 
distribution with mean equal to 0.5. The response Ui was generated from Normal(xii0ji -f 
Xi 20 i 2 +XiiXi 20 i 3 , 1), where On, 9i2 and 0^3 denote the true values used to simulate the data. 
We generated two clusters of observations of equal size, Si and 5*2, with ni = n 2 = 100, 
by setting: 0^.^^ = (3, 5, 9) in cluster 1 and 0si=2 = (O 5 ^5 0) cluster 2. 

The data generating process contained an interaction term only in one of the clusters 
{9i3 = 9). When htting the RPMS, SSM and the PSBP-MM we intentionally did not 
specify interaction terms in the regression sampling model. However the ability to perform 
variable selection jointly with covariate dependent clustering enabled the RPMS, the 
PSBP-MM and the PR to achieve robust predictive inference. To illustrate this property, 
let us consider the posterior distribution of the regression coefficients obtained by the 
SSM, RPMS and the PSBP-MM respectively. Given that a priori the cluster allocation 
of the SSM depends only on the cardinality of the clusters, the posterior of the regression 
coefficients under this model is invariant with respect to the different patterns in the 
covariate vector. Figure 1 presents the posterior density for the regression coefficients of 
the two covariates under the SSM. In the RPMS, since cluster allocation depends also 
on patterns in the covariate space, the distribution of the regression coefficients varies 
across different combinations of covariates. In our example there can be four different 
combinations. The fact that in RPMS the cluster assignment, and consequently the 
posterior distribution of the coefficients, depends on the covariates allows us to detect 
the effect due to the interaction term by inferring a cluster in which it is more likely to 
hnd both the covariates equal to one and then estimating the cluster-specihc regression 
parameters. This can be seen in Figure 1 (bottom), which shows the posterior density of 
the regression coefficients given that both covariates are activated. On the other hand, 
the SSM accounts for the interaction by estimating an extra component in the mixture 
distribution dehned for the regression coefficients (see top-left density in Figure 1). 

Similarly to what happens for the RPMS, the PSBP-MM assigns high probability to 
a mixture component which contains the combination of the covariates activating the 
interaction term. 

Obviously, this difference has a direct effect on the predictive distribution of the 
response. In Figure 2 we display the predictive densities for the four combinations of the 
covariates obtained when htting the SSM, RPMS, PSBP-MM and PR. 

The effect of the model misspecihcation becomes evident when looking at the predic¬ 
tive distribution for xi = 1 and ^2 = 0. It is worth noticing that, although the PR does 
not include a linear regression in the mean of the sampling model, this leads to robust 
predictive inference thanks to the model on the covariates. The latter permits to identify 
the four combinations of the covariates and then associates to each of them a cluster 
specihc mean in the response submodel. Consequently, PR identihes also the particular 
combinations of covariates that activates the interaction effect in one cluster. 
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Figure 1: Posterior density of 61 (left) and of 62 (right) in scenario 1 for SSM (top) and 
for RPMS (bottom). For RPMS, we consider the combination xi = X 2 = 1- 


PSBP-MM allows to select the covariates accounting for both the relevance in explain¬ 
ing the outcome and in partitioning the observations in clusters. We can summarize the 
importance of the d-th covariate by computing the quantity 1 — Pr( 7 id = ... = ''yKd = 
0 I y), i.e. the probability of inclusion of the d-th covariate in the model. The latter 
quantity takes a value very close to 1 for both covariates. This result is necessary for 
the PSBP-MM to achieve predictions robust to model misspecification, because it allows 
to capture the pattern in the covariates activating the interaction term in the sampling 
model. 

In addition, it is worth mentioning that grouping observations in clusters characterized 
by similar covariates may lead to identifiability problems of the regression coefficients 
within the model of the response in some clusters. The prior distribution over these 
regression coefficients together with the hyperprior distribution over the precisions of 
these priors allows very often to achieve robust predictive distribution. 

Figure 3 displays the posterior density of the continuous indicators employed by PR 
for performing variable selection. These highlight that both covariates are important in 
terms of determining the clustering structure. This is because the PR identifies clusters 
of response values sharing the same mean and the same combination of covariates, com¬ 
pensating in this way the model misspecification (note that in PR we are not regressing 
the response vector on the covariate matrix). 
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Figure 2: Predictive density of y for the four possible combinations xi = X 2 = 0, xi = 
X 2 = 1, (ail, 5 : 2 ) = (1,0) and (a;i,a; 2 ) = (0,1) in scenario 1 obtained fitting the SSM, 
RPMS, PSBP-MM and PR. The solid red line indicates the true density of the response 
for the four covariate combinations. 
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Figure 3; Posterior density of the continuous indicators tti and 7r2 for PR in scenario 1. 
Values close to 1 indicate the importance of the covariates for the clustering. 
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We conclude highlighting that we have not included an intercept neither for RPMS nor 
for SSM and the PSBP-MM and, accordingly, we have generated the observations from 
a regression model that does not include the intercept. This has been done to facilitate 
the presentation of the results. We have also performed the same simulation of scenario 1 
adding the intercept to these models and using observations generated from a regression 
model including cluster-specihc intercepts and we have obtained similar conclusions. 

6 Example: determinants of glycohemoglobin levels 

In this section we illustrate some of the discussed methods on a real data application 
aimed to identify the most relevant biomarkers of glycohemoglobin levels in diabetic and 
non-diabetic patients. Glycohemoglobin is hemoglobin combined with glucose and a high 
level of glycohemoglobin is usually associated to diabetes mellitus. Glycohemoglobin 
levels are measured as percentage of hemoglobin. Measurements of glycohemoglobin 
provide information on glucose levels over a period of three months, since once glucose 
combines with hemoglobin it can be traced for a period equal to the lifespan of red 
blood cells. As such glycohemoglobin is considered a better indicator of diabetes mellitus 
than direct measurements of glucose. Levels of glycohemoglobin higher than 7% are 
associated with diabetes, while average levels are between 4% and 6% in healthy people. 
Furthermore, high levels of glycohemoglobin are correlated with the risk of developing a 
variety of diseases such as diabetic nephropathy, neuropathy, angiopathy and retinopathy. 

The present section shows results obtained on a dataset containing 5089 patients for 
which the values of glycohemoglobin and of 22 covariates (9 binary covariates and 13 
continuous covariates) have been recorded. The dataset is available at www.biostat.mc. 
Vanderbilt .edu, and a description of the covariates is contained in Table 1. Incomplete 
records have been removed and income values have been discretized to three categories 
with cut-offs equal to $25000 and $75000. 

We compare the performance of Prohle Regression (PR), Random Partition Model 
with covariates Selection (RPMS) and Probit Stick Breaking Mixture Model (PSBP-MM) 
on this dataset and we focus mainly on describing the variable selection and clustering 
output. We have extended the PR and RPMS to include continuous covariates, by 
assuming Normal distributions within each cluster. A further modihcation of the PR has 
to be employed in order to perform variable selection on these continuous covariates, as 
it has been described in Section 4. 

6.1 Summarizing variable selection output 

The three models under analysis select important variables using different criteria. Ac¬ 
cording to PR the important covariates are those that contain clustering information. 
Liverani et al. [2015] propose to summarize the PR variable selection outcome through 
the distributions of vr^ | ^ (see Equation (19)), which has support on (0,1) and values 
close to 1 indicate that the d—th covariate is important. Figure 4, presents the posterior 
distributions of 7id \ y for all d and a quite large number of covariates seem important 
(posterior median of tt^ higher than 0.5), in particular covariates 1, 7 and 22 have a 
posterior mean between 0.5 and 0.7, covariates 4, 6 and 20 have posterior median be¬ 
tween 0.7 and 0.9 and covariates 2, 5, 8 and 9 have posterior median for vr^ larger than 
0.9. PR selects a large number of covariates because covariates can affect the value of 
the response exclusively through the clustering assignment. So, as it has been shown in 
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Table 1: Description of the covariates for the glycohemoglobin example. 


Number 

Description 

Unit 

Type 

1 

Income ($25000; $75000] 


Binary 

2 

Income >$75000 


Binary 

3 

Gender (male=l) 


Binary 

4 

Other hispanic (yes=l) 


Binary 

5 

Non hispanic white (yes=l) 


Binary 

6 

Non hispanic black (yes=l) 


Binary 

7 

Other race (yes=l) 


Binary 

8 

On insulin or diabetes medicines (yes=l) 


Binary 

9 

Diagnosed with diabetes mellius (yes=l) 


Binary 

10 

Age 

years 

Continuous 

11 

Weight 

cm 

Continuous 

12 

Standing height 

cm 

Continuous 

13 

Body mass index 

kg/m 

Continuous 

14 

Upper leg length 

cm 

Continuous 

15 

Upper arm length 

cm 

Continuous 

16 

Arm circumference 

cm 

Continuous 

17 

Waist circumference 

cm 

Continuous 

18 

Triceps skin fold 

mm 

Continuous 

19 

Sub-scapular skin fold 

mm 

Continuous 

20 

Albumin 

g/dL 

Continuous 

21 

Blood urea nitrogen 

mg/dL 

Continuous 

22 

Creatinine 

mg/dL 

Continuous 
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Figure 4: Posterior density of the continuous indicators tti, ..., 1122 in PR for the analysis 
of glycohemoglobin. Values close to 1 indicate the importance of the covariates for the 
clustering. 


the simulation study (Section 5), if the true relationship between the response and the 
covariates is linear, PR approximates it by dividing the covariate and response spaces in 
such a way to have in each part homogeneous values of both covariates and response. 

A different concept of variable selection is implied by RPMS, which performs cluster- 
wise regression and selects important covariates just for the linear model specified within 
each cluster. In this case we could summarize the global importance of the d—th covariates 
by the information contained in | ?/ (see Equation (18)). However, this information is 
quite hard to interpret because even for very high values of the d—th covariate can still 
affect the response by informing the clustering structure as well. For this reason Barcella 
et al. [2015] propose a two-steps approach which consists of first finding a posterior 
estimate of the partition of the observations, and then, conditional on such estimate, 
determining the posterior distributions of cluster-specific regression coefficients. Fixing 
the partition allows us to check the importance of the covariates in different clusters by 
computing the marginal probability of inclusion of each covariate within each cluster (see 
Figure 5). The covariates selected in the majority of the clusters are 8 and 9 followed by 
10 (which was not selected by PR), 20 and 22 . 

As noted previously, PSBP-MM does not assume any model on the covariates (see 
Equation (21)). Chung and Dunson [2009] propose a global null hypothesis for selecting 
the most important covariates. In particular, the d—th covariate is not important for 
the model if 71 ^ = 72 ^ = ... = 0. However, as acknowledged by the authors, such 
hypothesis is oversensitive given that the sequence of weights decays toward zero quickly 
and a covariate may start to be important just for very small mixture weights. So they 
propose to approximate the nonparametric model with a parametric version obtained 
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Figure 5: Posterior probability of inclusion of the covariates in RPMS for the analysis of 
glycohemoglobin. 


by truncating in Equation (21) up to some level K. We fix K=15. We evaluate the 
posterior probability 1 — Pr( 7 id = ... = 'yxd = 0 | y) for all covariates. Even if the 
results of RPMS for the covariates included in the majority of the clusters agree with 
those under PSBP-MM (with the exception of age), the latter selects also covariates 14, 
15, and 16. 

6.2 Summarizing clustering output 

The differences among the models outlined in the previous section affect also the clustering 
output. As models assume a random number of clusters, we first consider the mode of 
the posterior distribution of the number of clusters under the three models finding similar 
results for the models (10 for PR and 9 for RPMS and PSBP-MM). 

In order to get some understanding about the cluster configuration, we take as point 
estimate of pn \ y the configuration which minimizes the Binder loss function. PR cluster¬ 
ing is driven by different combinations of binary covariates, while all continuous covariates 
show similar patterns across clusters (except albumin and creatinine). Although RPMS 
includes the clustering information contained in the covariates similarly to PR, the clus¬ 
ters composition under the RPMS seems to be influenced equally by the discrete and 
continuous covariates. This is the result of having specified a linear regression model 
within each cluster which already accounts for the relationship between levels of the re¬ 
sponse variable and different combinations of the binary covariates. Finally, PSBP-MM 
does not account directly for possible patterns in the covariates and the clustering is 
exclusively in terms of the patterns in the response variable. However, the probability 
of the partition of the observations varies smoothly across the covariate space. This al¬ 
lows PSBP-MM to account implicitly for the effects of interactions among covariates, as 
shown in the first simulation study, by simply adding clusters. Similarly to RPMS, the 
compositions of the clusters seem to contain different combinations of binary covariates 
and continuous covariates. 


7 Discussion 

In this paper we have reviewed the most relevant literature to date on Dirichlet Process 
Mixture models with covariate dependent weights and the corresponding techniques for 
variables selection. Covariates can offer extra information on the partition of the data 
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and accounting for possible structure within the covariates can improve the predictive 
performance of the model. 

The most common solution to account for the presence of covariates within a DPM 
framework consists in assuming the covariates are generated by a probability distribu¬ 
tion and therefore specifying a joint DPM model on the augmented space including 
both response and explanatory variables. This solution is convenient in applications, 
as it becomes straightforward to obtain covariate dependent weights in the DPM model. 
Moreover, it leads to improved predictions in regression setups, still allowing for efficient 
computations. The major drawback of this approach is related to the fact that consider¬ 
ing a joint probability model for response and covariates may lead to the likelihood being 
dominated by the covariate specific terms, a problem that becomes non-ignorable when 
dealing with a large number of covariates. A solution to this problem is represented by 
the EDPM model (Wade et ah [2014]). This model introduces two clustering steps: first 
the observations are clustered on the basis of the response values and subsequently the 
model accounts for the covariate patterns within each of the clusters of the response. 

Alternative solutions can be found in the literature on covariate dependent Random 
Partition Models, in particular in research concerning the Dependent Dirichlet Process 
and dependent stick-breaking process in general. These techniques offer an elegant way to 
account for dependence of the weights in the stick-breaking representation on covariate 
information. However, when dealing with continuous covariates (or categorical with a 
large number of levels), these methods specify sequences of probabilities of cluster assign¬ 
ment for each observed level of the covariates, leading to difficulties in the interpretation 
of the clustering output. Moreover, posterior computations are often challenging when 
Gx is not marginally a DP, forcing the user to employ parametric approximations or 
expensive algorithms. 

The variable selection techniques proposed in the literature for covariate dependent 
RPMs aim at identifying the most influential covariate for the partition of the observa¬ 
tions. Especially for the augmented response models (and consequently also for PR), the 
likelihood of the model of the covariates can dominate the DPM when a large number 
of covariates is involved. By introducing latent variables we can eliminate the effect of 
specific covariates in determining the partition and consequently mitigate this problem. 

In many applications it is of interest to identify those covariates that best explain a 
response variable. The RPMS models extend the augmented response models to allow 
for variable selection in regression settings by specifying spike and slab distributions 
as base measures. Spike and slab priors are commonly used in the Bayesian paradigm 
to perform variable selection and recently they have been employed in non-parametric 
settings in context of DPM of regressions. We have shown through simulations that 
specifying a model on the covariates leads to inference robust to misspecification in the 
sampling distribution of the response. Of course, this comes at a computational cost. 
We have compared the performance of the RPMS with the SSM, a similar model where 
the covariates are considered fixed and not random. We have also presented results 
achieved using the PR and PSBP-MM. In the RPMS cluster assignment depends also 
on covariate information, while in the SSM it is affected only by the response values. 
This difference is reffected in the posterior distribution of the regression coefficients, 
which is dependent on the particular covariate pattern when fitting a RPMS. Obviously, 
predictive inference under the RPMS is conditional to the particular vector of covariates 
of a hypothetical new individual, while under the SSM the predictive distribution for 
a new individual is independent of his/her covariate profile. Also PR and PSBP-MM 
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are robust to model misspecification, thanks to the inclusion of the covariates in the 
distribution of the partition. Finally, we have also highlighted the different concepts of 
variable selection employed by PR, RPMS (or SSM equivalently) and PSBP-MM using 
a real example. 
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